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Abstract 

We investigate the behavior of a Bose-Einstein Condensate (BEC) under the influence of a central 
barrier as the particle number trends towards the thermodynamic limit. In order to perform these 
studies, we present a novel method which is tractable in the large- N limit. This method employs 
what may be considered to be a generalized Wannier basis, which successfully incorporates features 
of previous theoretical and computational assays to the splitting problem, including mean field 
effects, and has access to the dimensionality, trap parameters, and particle numbers relevant to 
recent experiments. At any barrier height we are able to discern between a two-mode state and a 
state which is described sufficiently by mean field theory and, further, give a criterion and technique 
for matching the two-mode theory to the zero-barrier state. We compare the basis used in this 
model to the de-localized basis functions underlying alternate models used in recent theoretical work 
on the double- well splitting problem and show that only the generalized Wannier basis displays the 
level crossing and emergence of two complex order parameters with overall U(1)(BU(1) symmetry as 
expected from a large- N analogue of the Superfluid to Mott insulator transition. Using this model, 
we identify a universal structure, independent of N, in this phase transition. We also present an 
analytic and model-independent description of this universal structure and discuss its consequences 
for realizing true two- mode physics with a BEC which trends towards the thermodynamic limit. 
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I. INTRODUCTION 

The realization of Bose Einstein condensation in a dilute gas of atoms [lj and the verifica- 
tion that the condensate order parameter both exists, and is characterized by the existence 
of a well-defined global phase, as witnessed by the observations of solitons j^j], vortices jsjjj^] 
and laser-like interference 5] has stimulated a great deal of research in recent years. Specifi- 
cally, a BEC coherently split in a double-well potential holds promise as a basic tool to study 
symmetry breaking, decoherence and phase diffusion properties of quantum systems as well 
as promise for use as an interferometric tool with an efficiency below the shot-noise limit jfj]. 
As such, these systems have attracted a lot of attention, with Saba et al. creating the first 

Q 

such BEC interferometer in 2004 [7]. These experiments have been refined with atom-chip 
technology to the point where a deterministic precession of relative phase js] attributed to 
differences in mean particle number in each well has been observed as well as a loss of a de- 
terministic interference pattern attributed to phase diffusion [9]. At smaller particle number 
Oberthaler et al. obtained fine enough control over a central barrier to reach the Josephson 

n 

regime in 2006 1_0J. 

One basic scientific question associated with the "splitting" of a BEC, first raised in Q| 
is whether, or under what conditions, the two moieties generated after splitting have the 
same phase and how long such phase-coherence lasts. Since a BEC is a mesoscopic system 
it is not a priori clear whether a classical or a quantum mechanical description of the order 
parameter is appropriate. In the former case, two independent condensates split from a com- 
mon progenitor would share a phase until interaction with the environment destroyed phase 
coherence similar to the way that the pieces of a cleaved crystal would retain a common 
orientation. In the latter case, a quantum description of the splitting process suggests, since 
relative number and relative phase are conjugate variables, phase coherence will be destroyed 
for independent condensates. Lattice experiments operating at small particle per lattice site 
exhibit a loss of phase coherence during the so-called Superfluid to Mott insulator tran- 
sition, but show a restoration of phase-coherence faster than so-called "phase-incoherent" 
states once the lattice depth is decreased [llj. It is not clear, without a theoretical tool 
to accompany these experiments, what is responsible for these two phenomena, where both 
the Mott Insulator state and the phase-coherent state are described as "fragmented," or 
pure Fock states with perfectly defined particle number. Such a theoretical tool should 
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also be able to investigate what barrier heights and ramping times are needed to engineer 
"squeezed" and other exotic states. We introduce such a theoretical method capable of ac- 
cessing physical data inaccessible to in situ imaging in order to compliment the set of BEC 
splitting experiments currently being performed in hopes of answering these fundamental 
and technical questions. 

Theoretical descriptions of certain aspects of degenerate bosonic systems in one or more 
modes have been contributed from a variety of disciplines within physics over the last decade. 
The description of a BEC in terms of an order parameter with definite phase as an example 
of spontaneous symmetry-breaking is inherited from the condensed matter literature, as is 
the use of model many-body wavefunctions and Fock-space expansion coefficients to describe 
particle distributions in between multiple wells, aka the Bose-Hubbard model. When de- 
scribing a two-mode or splitting BEC, Bose-Hubbard type models use tunneling rates and 
site energies, which are typically input as parameters or approximate wavefunctions [12j but 
may actually have a complicated dependence on underlying many-body wavefunctions and 
trap geometry as well as having dynamics of their own during the splitting process. More 
realistically, in the context of a highly-restricted "quantum optics" type quasi- Gaussian 
ansatz, Zoller et al. [13| investigated general properties of the dynamical splitting of a BEC. 
Perhaps most importantly, they noted that because there are two different types of dynam- 
ical variables of interest in the system: the spatial variables which govern how the density 
of the atomic cloud(s) evolve and the Fock space variables which govern whether particles 
are localize into a single well, there are two different types of adiabaticity, each associated 
with a characteristic timescale. We further note that since in situ imaging of a BEC only 
probes the atomic density of the cloud, the distribution of the Fock space variables of the 
system must therefore be inferred from the density, for instance by the loss and revival of 
an interference pattern. The detailed analysis of one number-squeezing lattice experiment 
[jjj], however, has shown that mean- field effects may mimic the experimental signature of 
number- squeezing [ijj] indicating that in situ imaging should not be used as a stand-alone 
diagnostic of the Fock-space distribution of a split or squeezed BEC. However, an ab initio 
computational method allows this variable, the primary component in engineering quantum 
degenerate states beyond mean-field theory, to be investigated directly, b y c onstruction. 

Recently, more exact, "multi-configurational" schemes If]] 17] is| 19 1 20], derived from 



methods familiar to quantum chemists and many-body nuclear theorists have been developed 
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which self-consistently include both mean-field and particle number dynamics for a set of 
indistinguishable bosons with access to more than one mode. In spite of all of these efforts, a 
single theory which synthesizes the relevant aspects of these approaches and calculates useful 
quantities like phase diffusion rates and the effect of barrier raising times on the subsequent 
interference patterns has not yet been presented. 

There are two essential difficulties in correctly simulating the splitting of an initially- 
coherent BEC into two distinct moieties. First, is the fact that as the barrier raises, the 
system acquires an additional degree of freedom which was absent at t = 0. Namely, at 
some point in the splitting, a breadth of Fock space expansion coefficients corresponding 

, is needed. Conversely, mean field 



2 if ] 22] - which is the correct theory 



to the possible distributions of particles in either wel 
theory, as described by the Gross-Pitaevskii equation 
at low barrier heights - assumes that all particles are contained in a single quantum state. A 
successful theory, therefore must correctly discriminate in between a one-mode (GP like) and 
two-mode states in the course of a time-dependent and potentially non-adiabatic process. 
The method described in this paper applies the Penrose-Onsager criterion to a form of the 
reduced density operator following [13J in order to do this. The second difficulty is that, if 
one takes a condensed matter theorist's viewpoint, the complete fragmentation of a BEC 
into two independent BECs should also be an example of symmetry breaking, in which 
a system described by a single order parameter with a global U(l) symmetry, becomes 
a system with U(l) © U(l) symmetry. That is, each of the two wells should contain an 
independent condensate such that an arbitrary phase shift on one should not affect the total 
system energy. In order to motivate our choice of spatial basis functions, we show that this 
latter requirement is not satisfied by a naive implementation of even the sophisticated self- 
consistent methods such as in [23]. In response to these difficulties, we present a basis and 
method which smoothly transitions from a single totally-occupied zero-temperature ground 
state to a fragmented state in such a way that the initial density as well as the final densities 
and Fock-space distributions have physical meaning. 

In Section [Til we develop such a theory which is correct in both the low and high-barrier 
limits. First, we develop equations of motion for two spatial mode functions and the possible 
partitions of particles in between them which are tractable and separable in the limit that 
N ^> 1. We give a criterion for when the system is described by a single totally-occupied state 
despite using an explicitly two- mode basis and give a further means of checking/correcting 
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all necessary physical data of the system in that case. In Section HTT1 we briefly describe the 
numerical implementation of the method described in [Til 

In Section IIVI we show results from a set of calculations which verify that the basis 
we have employed, the generalized Wannier basis, correctly reproduces the Superfluid to 
Mott insulator transition (SF-MI). In order to do this, we define what it means within 
this theory to be in a Superfluid state and within a Mott insulator state, then we show 
that a curve crossing exists in which the energy of the Mott insulator state eventually 
falls below that of the Superfluid state as a function of barrier height. Conversely, we 
show that in the same physical system treated with a de-localized basis, the Superfluid state 
always lies energetically below what is defined as the "fragmented" state. Therefore, no such 
curve crossing exists within the analogous "single-" and "double-macroscopically" occupied 
gerade/ ungerade states. Finally, in order to understand how these specific computations 
will generalize to other systems, we present estimates of how the band-gap of a de-localized 
basis scales with particle number, trap frequency and interaction strength. We show that 
this is strictly positive for repulsive condensates and, as such, the SF-MI transition cannot be 
simulated in models employing such a de-localized basis. On the other hand, an analogous 
calculation for the generalized Wannier basis model shows that, at high barriers, a Mott 
insulator phase is always the ground state for repulsive interactions. Therefore, the curve 
crossing, and hence the SF-MI transition, is a universal feature of this model. 

Our initial physical results from this theory, described in |Vl follow from a set of numeri- 
cal and analytical investigations of the splitting process as the number of trapped particles 
trends towards the large- N limit. We are able to delineate the regime in which a two- mode 
model is appropriate and this analysis indicates that there is a very narrow region, charac- 
terized by a universal mathematical structure, in which two-mode models are applicable to 
splitting process. 

In Appendix A a method to generate an effective ID equation which self-consistently 
incorporates trap and mean-field data from the transverse directions when there is only one 
principal (splitting) axis of interest is described. 
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II. THEORY 



Here we develop the basic dynamical laws which we use to describe a set of N identical 
bosons in an external trapping potential V ext (r) and give the means to interpret the state 
vector. 



A. State Vector and Equations of Motion 

Starting from the second-quantized Hamiltonian in the contact-approximation for quasi- 
1D (see |A] for a discussion of how we self-consistently include transverse trap data when only 
one splitting axis, r, is of interest): 

H = Jdr[¥(r)(f + V ext {r))i>(r) + ^¥(r)¥(r)i>(r)i>(r)l (1) 

where g gives the interaction strength through the s-wave scattering length, a s , as g = 
47r ^ r ' 2 as is appropriate for a dilute low-temperature gas and the & are the bosonic field 
operators satisfying the usual commutation relations [^(r),^(r')] = [\&t(r), ^t(r')] = 0, 
and [^(r), ^(r')] = 5(\r — r'|). Because we are including the possibility that the state is 
a correlated two- mode or fragmented state, we must use two-mode field operators \l/(r) = 
ai0i(r) + a 2 </> 2 (r). 

To describe all possibilities for the occupation of these modes, the full state vector is 

n v a\ \/{N — a)\ „ 
At this point the standard two-mode operator algebra for bosons [12] gives 

2 2 

E = ($\H\$) = ^ Pi k€ i k + f Pjklm^jklm, (3) 

j,k=l j,k,l,m=l 

with the auxiliary definitions p^ = ($|aja fc |$), pjkim = ( < ^ ) |Sja|.a;a m |$), 

e jk = J rfr ^K~^~~~ v2 + V ext)4>k: and T jHm = J dr<f>j(r)<j>l(r)<j>i(r)(j) m (r). 

We now construct the Action in order to generate equations of motion for all of the 
quantities declared to be dynamical variables 

S = f dt{($\H -ih§T-J2 ^ k I dr ^ k - ( 4 ) 
J j,k=i J 



where the Lagrange multipliers fijk are introduced to enforce the constraint J dr<p*(j) k = 5 
Using the above 



S= f dt{E - ih{$\^-\$) - fa* I drfifa - 5 jk )} 
J j,k=i J 



(5) 



we now declare the set of Fock-space coefficients {C a } and the spatial functions {4>j{r)} 
to be dynamical variables. Employing the Dirac-Frenkel variational principle, as is done for 
the case of arbitrary modes 18J, the equations of motion for the Fock-space coefficients are 



given by the condition -j^ = gives 



dS dE 

— ^- - % 



N 



dC, 



dC* 



(6) 



2 2 
j,k=l 7 



2 S £U dCt < 



jklm 



or, in a more compact notation 



N 



0=0 



where 



Kip = (7l { X] e Aa]a k ) + | X r iWt»( ft ''' ' 



(7) 



(8) 



/ Q fcJ + n L jkim{a'ja k aia m ) f | >,). (9) 

>j',fe=l j,k,l,m=l 

Similarly, the equations of motion for the two mode functions are given by the condition 



OS 



i*(r') 



0. Now, using the fact that Q(j) i^ = ^j q S(\r — r '\)i this condition yields: 



j,k=i 



j,k=l 



Pjklm 

j,k,l,m=l 



v f.->,.OlO„, + (pjOkqOiU,,, I 



^2(fX jk S jg (l) k ) (10) 

j,fc=l 



ifr^Pqktfik = Yp q k(f + V ext )(f) k + g ^ Pqklm(<Pt < Pl < Pm) ~ /XVqk&k) ■ (H) 
fc=l fc=l fc,i,m=l fc=l 

These variational calculations have been performed elsewhere in the context of identical 
and distinguishable bosons [181] [24| up to this point. Typically, the derivation continues 
by multiplying either side of the above equation by (pqk)^ 1 in order to decouple the two 
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time derivatives, however, this produces singular equations of motion when p q k has a zero 
eigenvalue (i.e. when only one state is occupied). Since this is precisely the initial state 
we wish to consider in the case of a single well deformed into two wells, another procedure 
is employed. We will observe and comment on some ramifications of using schemes which 
become singular in the case of only one occupied mode in Section IIVI 

In deciding how to proceed we need to know what approximations to the state vector are 
relevant in the large N limit. We do not parametrize, truncate or approximate the Fock- 
space variables keeping the entire expansion (J2]) and the exact equation of motion (JSJ), since 
the degree of coherence and fundamental interpretation of the system will depend vitally on 
both the breadth and relative phases of the distribution of the variables {C a }. 

In order to decouple the equations of motion for the {4>j}, we note that the approximate 
form of the distribution of the {C a } is binomial for the case of any ground state configuration 
of the symmetric double- well if the mode corresponding to <p\ (02 ) is approximately left 
(right) localized. 

More strictly, in the case of repulsive interactions (g > 0), the following bound holds: 



1 N/2 / AH 

lCal< l v SRi' (12) 

Therefore, the standard deviation of the distribution of the {C a } should scale as ~ ^= and 
the numerical prefactors on fill I) which derive from the quantities pjk, Pjkim will be well 
characterized by their Fock-state values when N ^> 1. Consequently, in the thermodynamic 
limit, the equations of motion for the lobes are given, to leading order in 1/N, by the 
"diagonal" contributions: pjj, Pjkjk, Pjkkj- 

Using this large- iV approximation, ( II ip becomes, explicitly writing out the components: 

Pix{f + V ex t) 

P22(f + V ext 

I Pllll|0l| 2 +Pl22l|02| 2 P12120201 , . j (10) 

Y P21210102 P2112|0l| 2 + P2222I02I 2 

In the above, we also use the fact that, for unitary time evolution, the norm-preserving 

Lagrange multipliers, pn and p 2 2, are unnecessary. 

As the art of variational science is to find a variational space that gives the desired 

limiting cases automatically, we show that, by using as <pj what may be considered to be 
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generalized Wannier functions satisfying the non-linear equations (|T3l) . we get the correct 
high barrier limit of two fragmented condensates. We discuss this in detail in Section IIVI 
giving a comparison to other assays at the splitting problem. 

Finally, we point out an external criterion with which we can ensure that our theory is 
consistent with correct theory in the low barrier limit (i.e. mean field theory as encompassed 
in the condensate order parameter and Gross-Pitaevskii equation). For a zero-temperature 
Bose-Einstein condensate, all physical data resides in the complex order parameter $op(r, t) 
which satisfies the relationship §op( r ,t) = \/ Ps( r , £)e ie(r,i ), where ps(r,t) is the spatial 
density of bosons 25[. As such, even in the case when a compound or parametrized com- 
putational basis is used, the above relation can be used to check the consistency of our 
ansatz. 

We further use the sufficiency of the complex order parameter to determine a renormalized 
interaction strength as a means to compensate for any errors introduced by the possibly 
"spurious" kinetic energy terms in the generalized Wannier basis or from using the uncoupled 
Fock-state equations for the {4>k}- The procedure we adopt is to use a diagonal and off- 
diagonal interaction strength, g and gu respectively, in order to describe the strength of 
the particle-particle interactions. To keep consistency with the high-barrier limit of two 
independent condensates, the |0fc| 2 0fc terms use the prefactor g, however, we mediate the 
strength of the \(f)j\ 2 (f)k (k 7^ j) by a g 12 chosen to give the correct density, and therefore 
order parameter, at zero barrier. While this procedure introduces an ad hoc modification to 
the equations of motion, we note that, because of the sufficiency of the order parameter, this 
does generate the exact ground state of the system when correctly interpreted. In practice, 
we find that the modulation of gu is on the order of a few to tens of percent and has no 
significant effect on the basic scientific conclusions derived from this computational theory. 

With a final rearrangement of terms and introduction of the renormalized parameter 
#12, Our final equations of motion, a coupled system of equations for the Fock and spatial 
dynamical variables, therefore, appear as: 

N 

ihc^ = J2n^Cp (14) 

/3=0 

for the {C a }, with % 1 p defined in ([HI) and 
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' 01 j 







/il2 ' 


II 




A«2i J 




\<h ) 



Pn(T + V ext ) 

P22(f + V ext ) 

+ | fi'l0l| 2 Pllll 25-120201P1212 , , , . 
I 2#i20i02P2121 # 1 02 1 2 P2222 

Where in the above the renormalized parameter gu has been introduced and the identity 
P1221 = P1212 has been used to allow a useful rearrangement of terms. 




B. Interpretation of State Vector 

Despite explicitly using two spatial basis functions {01,02} in this theory and a full 
breadth of Fock-space variables {C a }, we are interested in determining whether the system 



is characterized as two-mode vs. one-mode. In order to do this, we follow [13( and use the 
Onsager-Penrose criterion. By tracing over the spatial variables of the density operator one 
gets a 2 x 2 matrix. 

PF = { {aW) (a[a2) ] (16) 

This reduced density operator has two eigenvalues A + /_ which sum to N. In the case that 
A + = iV and A_ = 0, the system could be sufficiently described by the GP equation. 
Should the Fock-space density operator yield eigenvalues A + = A_ = N/2, then the system 
is "fragmented" into two independent condensates. Finally, a correlated two-mode model, 



such as has been used to describe Josephson junctions 



261 ] 271 | and BEC atom interferometers 



28] 29], is only appropriate when > A + > N/2. 

The other reduced density matrix, i.e. the density matrix traced over the Fock-space 
variables, returns the spatial density of the system and is given explicitly by 

Ps = {aW)^ 2 + 2^{{a\a 2 )4>\4> 2 } + (a 2 a 2 )|0 2 | 2 (17) 

this is the quantity we use, along with the Landau's identification, which becomes: 

$op(r,t) = V^M)e i0 ^ (18) 

in order to verify that our two-mode basis correctly reproduces all the physical data of the 
mean-field state when the mean-field description of the condensate is appropriate. 
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III. IMPLEMENTATION 



In this section we describe how we solve and implement the theory described above. 
Because this paper, the first of a two-part series, will be concerned with comparing the 
stationary states of various theories, this amounts to diagonalizing the equations of motion 
described above for the various basis sets under consideration in Section UVl Time evolution 
and applications of this formalism to barrier raising and ballistic expansion will be the 
subject of Part II. We find the ground state of the system by employing a relaxation method 
to the complex time version of the equations discussed above. This is done by taking the 
Wick-rotated equations of motion t — > r = —it and adding self-consistent estimates of 
the chemical potential and system energies from the spatial and Fock equations of motion 
respectively. The Wick-rotated equation of motion for either variable are: 




Pn(f + V ext ) \ ( 0i \ / /in /ii2 \ / 0i 

l>->-Al ■ \'<rl) I \ 02 / \ /<2J 1 1 22 / \ < >2 




h 



. g\<t>l\ Pllll 2#i2020lPl212 . . , . 
2#120 : I:02P2121 #|02| 2 P2222 

ar N 

h ^ = E °-J2 n ^ ( 2 °) 

/3=0 

Now, instead of the evolution characteristic of the Schrodinger equation (j>j(t) oc e~ ltfM / h 
and C a (t) oc e - ltE/h , the equations dUDO^n]) evolve form of exponential decay when the 
state has higher energy /chemical potential than the true ground state: (f>j(r) oc e^^" 1 ^^ 
and C a (r) oc e T ( E o- E )/ tl g that, as the quantities fijj and E iteratively converge to the 
ground state of the chemical potential and system energy, respectively. This is done by 
updating these quantities by the identifications fijj = ihjdx (j>*4>j and E = ^ C^H a ^Cp, 
for the provisional values of and C, renormalized after a few e-foldings of the relaxation 
process. 

When considering the excited states of the gerade, ungerade basis discussed in Section 
IIVt the configuration of C a for each such excited state is known a priori, and in this case, 
for a fixed distribution of C a , the system f fl9|) will unambiguously converge to the functions 
and chemical potential for that excited state and not the true ground state of the system. 
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One technical point is that since the system of equations comprised of (|T9|) . and (!20|) is 
non-linear, there is no guarantee that a stationary state produced by this method is a true 
ground state of the system. Stated another way, non-linearity means that there is no single 
energy surface for these variables since the energy surface depends on the {4>j} themselves. 
In order to ensure that we have not converged to such a "non-linear local minimum," we 
perform this relaxation method on several different guesses for the initial conditions and 
verify that they converge to the same answer. 



IV. BASIS COMPARISON 



A. Numerics 



The typical choice in discussing the splitting problem is to select basis functions of gerade 
and ungerade symmetries and use them to create approximately left and right-localized 
functions in order to describe the high barrier state. This was first done in the context of the 
statics of a weakly-interacting gas in a double well potential in [l2| using the ground (4> g ) and 
first excited state (<f) u ) solutions of the Schrodinger equation and their linear combinations 
0l = + and (pn = ^((fig — cfiu). This has the appeal that, in the non-interacting, static 
limit, at zero barrier + 4>r) n = 4>gi the true many-body ground state of the system. 

More recently, the dynamics of the splitting process has been investigated using the 
{<Pg-i 4>u} basis directly with the physical identification that if both gerade and ungerade 
functions are macroscopically occupied, the system is "fragmented" |23l||30]. 

We use neither of these methods and, in this section, point out some pathologies associ- 
ated with either the gerade, ungerade scheme. Instead we show that a basis which can be 
considered to be the non-linear generalization of Wannier basis functions generates a theory 
with a level crossing as the central barrier increases and show that this is a generic feature 
of that basis. 

We show that no such level crossing exists in the {<f) g , <p u } basis for the same system and 
give a criterion for when the ungerade states become important as the system trends towards 
the thermodynamic limit. 

Considering a system of 200 87 i?6 atoms confined in a trap with frequencies u x = 2n x 
44.7Hz, Uy = u z = Co>j_ = 2tt x 1.1kHz in the presence of a central Gaussian barrier 
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(Un)gerade basis functions and potential at various barrier heights 




-04 1 1 4 

10 15 20 25 30 10 15 20 25 30 10 15 20 25 

position (jim) position (jim) position (^m) 



Generalized Wannier basis functions and potential at equivalent barrier heights 




-0.2 I 1 1 1 1 -0.2 I ' 1 ' 1 -0.2 1 1 ' ' 

10 15 20 25 30 10 15 20 25 30 10 15 20 25 30 

position (urn) position (urn) position (^m) 

FIG. 1: Picture of basis functions used in [23|] and this work 

characterized by a variance of 10/zm, we compute the energies of various configurations of 
interest. The complex time versions of the effective ID equations of motion developed above 
in sections II I II and fX] are implemented in order to perform this study. 



B. Analytic model of energy crossings 

We present a discussion, in the context of this model, of generic features of either choice 
of basis and show that the Wannier basis always has a Mott Insulator ground state at high 
barriers while the gerade, ungerade does not. Further, it can be shown that the ungerade 
state becomes irrelevant in the large- N limit, as the energy gap scales as gN - an alternate 



proof of the perturbative result of Huang and Yang [31]. 

In order to show that the Mott Insulator state is a generic feature of the Wannier basis, we 
show that the Hamiltonian expressed in the Fock basis ([9]) has a minimum at N\ = N/2 and is 
concave. In the high-barrier limit, as in lower-right panel of FigHJ when the basis functions 
are localized / <ir |0i| 2 |0 2 | 2 — > 0, the (otherwise pentadiagonal) Fock-basis Hamiltonian, 



becomes diagonal: 



2 2 

= (01 { ' -a/ ( "/'./) + 1 \ ( 21 ) 

. 3=1 j=l 
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a) Lack of level crossing in g/u basis 5 b) SF-MI transition in Generalized Wannier basis 




w Barrier height (arbitrary units) Barrier height (arbitrary units) 

FIG. 2: (Non-)existence of curve crossings in de-localized and Wannier basis theories. In a) the 
heights of various excited states (as defined by promoting an additional particle from g to u) above 
the ground state (all particles in g) is shown. No crossing is observed in this basis, indicating that 
observed population of low-barrier "excited" states in g/u models may be a function of equations 
of motion which require initial population of the ground state. In b), we observe an energy level 
crossing consistent with the SF-MI transition in the Wannier basis. In this panel, the two-mode 
coherent state is defined by a Binomial distribution of \C a \ 2 and self-consistently relaxed and 
the "fragmented" state is a single, unit C a= N/2 as described in the text. 

and the energy becomes 

N 

(H) = c l uc ? = \ c p\ 2 (ciiP + e 2 2 (iV - p) + |r im /3(/3 - 1) + |r 2222 (iV -0)(N-P+ 1)) 

13=0 

(22) 

Invoking the symmetry of the well, r im = T 2222 for instance, the energy then can be written 

N 

{%) = £ C;HCp = \C p \ 2 (e n N + |r mi (2/3 2 - 2iV/3 + N 2 - AO) (23) 

0=0 

Which is concave and has a minimum with respect to /3 at = N/2. Finally, since H 
is diagonal when the functions <pj don't overlap, the ground state is a single configuration 
C/3 = 5pji/2, or Fock state and the system is in a Mott Insulator phase. 

Conversely, there is no analogue for this proof for the de-localized functions {4> g ,<j) u }, or 
their nonlinear generalizations. We can, however, examine how relevant the excited states, 
as defined by the promotion of an additional particle from cf) g to <p u are in the large iV 
limit. Looking at the energy difference between all particles in (f) g (since the quantities in 
(J9j) depend on particle number, we denote the ground state factors of e and V by o) and 
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N — 1 particles in <p g one particle in the first excited state (quantities denoted by +): 

AE = (H) + - (H) = (e+ + (N- 

+ |((iV - l)(4r+ g J + (N-1)(N- 2)T+ gg ) - (Ne° gg + 9 -N{N - l)T° gggg ) 

Since we are simply looking for a scaling argument we can equate the geometric quantities 
which only differ by the exclusion of a single particle (for instance ^ ggg = ^°ggg g ) i* 1 the limit 
of N 3> 1 giving. 

AE = e uu - e gg + |(4(iV - l)T gugu + (2 - 2N)T gggg ) (24) 
Roughly, one can estimate this by noting that T gugu and T gggg ~ 1 and e uu — e gg ~ fko giving 

AE ~ hw + ^A^ (25) 



A complimentary derivation of the energy gap for a dilute Bose gas first published in 31] 



and strong indication that the energy crossings and population of excited states reported in 
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30| may 



computational attempts on the splitting problem which use a de-localized basis 
only be relevant at low particle numbers and those results may not generalize to large- N 
systems. We further point out that using these methods to compute factors of T and e exact 
spectra for arbitrary trapping potentials and particle numbers can be calculated when the 
de-localized basis is used. 

We make two remarks here about the presence of this energy gap in models which employ 
a de-localized basis. First, this calculation shows that the alternate definition of "fragmenta- 
tion" which has appeared in the literature, namely, that a condensate is "fragmented" when 
both the g and u states are macroscopically occupied is not equivalent to the definition of 
a single unit-probability partition of the particles in between wells which is native to Bose- 
Hubbard type models. The former definition corresponds to a state which is, generically, 
of more energy than the ground state of that theory (all particles in g) and therefore must 
be reached by a non-adiabatic process. Conversely, we have shown that the Mott insulator 
ground state is a generic feature when Wannier basis methods are used to model repulsive 
condensates. 

Secondly, in the context of dynamical investigations of the splitting process, the non- 
invertability of pjf. (as discussed in Sec. [TTJ) means that in order for the equations of motion 
such as those used in [3] and 24] to be nonsingular, the initial state must include some 
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finite population of an energy level above this gap which is increasingly unphysical in the 
large- N limit. We note that if the structure of the energy levels seen in Fig. [2] a) is general 
and the energy gaps maintain approximatel y-eq ual spacing as the barrier raises, this may 



explain the oscillatory behavior observed in 



23|, for instance. 



As a graphical summary of the physical interpretation we propose for our theoretical 
description of the splitting process, we include figure [3j In this figure, we again find the 
ground state of the theory outlined in [Til and II III for the 200 87 Kb atom system described 
above. This figure shows relevant physical data for a range of barrier heights. Clearly 
visible in Fig |3k) are three distinct regimes, a coherent regime in which all N particles 
reside in a single state and GP theory should be applicable; a two-mode state, where there 
are two populated eigenvalues and it is sensible to talk about relative populations of two 
mode functions; and finally a regime in which there are two "fragmented" BECs which are 
described by who independent complex order parameters. 

As discussed in the figure caption, the eigenvalues show how to interpret the state. When 
there are {N, 0} eigenvalues, only the final density (and phase, should one be doing dynam- 
ics) has a clear physical interpretation through the Landau identification ( fl8i) and the C a 
and are simply elements of a computational basis. On the shoulder of the best-fit tank 
functions, a two-mode model is appropriate and the <pj become the mode functions. In 
this case, the C a also have a direct physical interpretation as \C a \ 2 being the probability of 
detecting a particles in <\>\. Finally, when the system has two N/2 eigenvalues, only one C a 
has non-zero probability and the system is correctly described as two independent conden- 
sates. If this were at the end of a dynamical splitting process, one would say that the initial 
condensate has been "fragmented" into two condensates. 

We have shown that the basis we have selected to generate the theory of a splitting 
BEC correctly displays a curve crossing in energies corresponding to the Superfluid to Mott 
insulator transition and that this is a generic feature of the model. We have also shown three 
distinct physical regimes of the splitting BEC and have carefully delineated what parts of 
the computational basis have a clear physical interpretation in each of these regimes. 
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(a) Eigenvals of Fock-space density operator p F 
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FIG. 3: A selection of the physical data which comprise our theory. In a) the system eigenvalues 

are shown as a function of barrier height, clearly spanning three distinct regimes. In b) the 

computational variables from zero-barrier state is shown, only the final density has direct physical 

meaning. In c) the data from the two-mode regime is shown, both the (f)j(r) and C a have the 

interpretation of the spatial and Fock-space variables of a two-mode theory. In d) the data for 

a "fragmented" state is shown in which both 4>i and (j>2 have the interpretation of independent 

complex order parameters and there is no relative particle number uncertainty. 
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V. ANALYTICS AND NUMERICS OF THE RELEVANCE OF A TWO-MODE 
MODEL 

We use the formalism developed in section [TT] to determine the ground state of a variable 
number of bosons in the trap discussed above. Our goal is to understand how various aspects 
of the splitting process scale with the total particle number N and what generic features can 
be observed. We observe that a type of universal behavior is associated with the splitting 
of a dilute Bose-condensed gas and we present a simple analytic model to understand this 
behavior. Finally, we discuss the ramifications of this universal feature of a splitting BEC. 

A. Numerics 

As in the preceding section, we use the complex-time form of the equations shown above 
and a relaxation algorithm in order to compute the ground state of a system of various 
values of N 87 Rb atoms confined in a external potential characterized by trap frequencies 
uj x = 2n x 44.7Hz, u y = u z = u± = 2n x 1.1kHz and a central Gaussian barrier, applied 
along the x direction (See |A] for details of how this system was modeled as an effective ID 
system) . 

In accordance with discussion in [III i n order to discern in between a single-mode "GP- 
like" state and a system that truly requires a two-mode description, the quantity of interest 
is the density operator traced over the spatial degrees of freedom or the "Fock-space density 
operator." ( !T6|) . We plot the eigenvalues of (|T6|) as a function of increasing barrier height for 
four decades of particles, N = {20, 200, 2000, 20000}, in the trapping potential described 
above. The results of these calculations are shown in FigHJ 

In the context of the splitting of a BEC as a phase transition, this figure may be under- 
stood as the left most region of each graph, characterized by eigenvalues of A + /_ = {N, 0} 
as being a "GP-like" state, namely, the Gross-Pitaevskii equation is sufficient to describe 
the zero-temperature physics of this state. At high barriers, when the state is characterized 
by A + /_ = N/2,N/2, the two generalized Wannier functions are separated enough that an 
arbitrary phase shift on one of them, will not affect the total system energy, satisfying the 
requirement that this process be a type of phase transition. The region in between these 
two asymptotes, or phases, is where it is appropriate to use a two-mode model to discuss 
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FIG. 4: Eigenvalues of the Fock-space density operator as a function of barrier height for four 
decades of N. In each case, each of the three regimes discussed in the text corresponds to a finite 
range of barrier heights. 

the relevant physics. 

Since the production of the first dilute alkali gas BEC, schemes have been proposed in 
which atom interferometry is performed on an initially phase-coherent two-mode BEC, see 
for instance 32J 33], therefore it is of fundamental interest to understand when a two- 



mode model is relevant. We devote the rest of this paper to understanding under what 
circumstances a single coherent, BEC can be brought into the two-mode regime. 

First, in order to examine the splitting process in a system-independent manner we define 
the intensive quantity we call the entanglement of the system which is simply equal to the 
normalized difference in between the eigenvalues of ( |T6l) : entanglement = (A + — A_)/./V. 

In Figure [5] we show the entanglement plotted versus the height of the central barrier 
for a few decades of total particle number, showing the transition from a unit- entanglement 
"GP-like" state through a system correctly described by a two-mode model and finally to two 
independent condensates described by entanglement zero. We further roughly scale out the 
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Entanglement parameter vs. Scaled barrier height 
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FIG. 5: Scaled graphs of the entanglement which suggest a lack of non-analyticity in the splitting 
process as N — > oo. 



dependence of the barrier height for different particle numbers by plotting the entanglement 
versus Vbarrier(% = L/2)/[j,n, where /jln is the chemical potential of each system, N = 
{20, 200, 2000, 20000}, at zero barrier. 

While it is typical to expect some non-analyticity to emerge in the thermodynamic limit, 
Fig suggests, however a universal feature associated with the process of transitioning from 
a one-mode state through a two-mode state in order to fragment a BEC. We conclude that 
this is, in fact, the case and present an analytic model to supplement this numeric data in 
the next subsection. 

Since the peak barrier height, even when scaled to the chemical potential, is not inde- 
pendent of the system geometry, we show how the entanglement of the N = 20, 000 particle 
system depends on the ratio of the self-energy to the hopping energy e^/en. These quanti- 
ties, typically denoted J and U, are analogous to the on-site and hopping strengths of lattice 
Bose- Hubbard type models [341] . Figure [6] shows that, while a two-mode model is appropriate 
over several decades of the hopping rate before fragmentation, a single-mode GP-like state is 
a necessary and sufficient description until the ratio of on-site to hopping strength, (e^/en) 
is well below 10~ 6 . In addition to the entanglement, we plot the expected phase uncertainty 



of a two-mode model using the Heisenberg relation AiVA0 1 35( as a rough guide to 



when two-mode interferometry would fail for lack of a well-defined relative phase. 
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Entanglement and Phase Uncertainty for N=20,000 
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FIG. 6: Logarithmic scale graph of the entanglement as a function of (unitless) tunneling rate for 
N = 20, 000. 

B. Analytic model of splitting regime 

In order to understand this behavior, we present an analysis of the basic mathematical 
quantity which dictates splitting process by assuming a simple analytical form for the dis- 
tribution of the {C a }. Examining the Fock-space density operator ffT6l) . for a symmetric 
geometry, the diagonal entries are known: (a\a\) = {a\a 2 ) = N/2. Therefore the entangle- 
ment is determined solely by the off-diagonal entry p 12 = (a\a 2 ) and its complex conjugate. 
Explicitly, this is quantity is given by 

N N 

p 12 = (a\C* a (a\a 2 )Cp\(3) = ^ C*C a ^y/N - a + lyfc (26) 

a,0=O a=l 

When this quantity vanishes, the state is fragmented or two independent condensates and 
when it equals N/2 the system is at unit entanglement, and a description by the GP equation 
is appropriate. The sum in (12 61) is a weighted, unit-offset autocorrelation of the distribution 
{C Q } and, by assuming an physically-motivated analytic form of {C a }, we are able to get 
an insight into the process of fragmenting a BEC. 

(a-N/2) 2 

Assuming a continuous, Gaussian probability distribution of \C a \ = ^=f e > 
the limits discussed above are represented by means on the variance of the distribution, 
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FIG. 7: Graph of p\2(2/N) vs. width of the Gaussian ansatz C a distribution spanning four decades 
of N showing the regimes of entanglement a) N = 20 87 Kb atoms, b) N = 20 87 Rb atoms, c) 
N = 2000 87 Rb atoms, d) N = 20, 000 87 Rb atoms. 
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FIG. 8: Universal feature of (2/N)pi2 = entanglement as desribed by the Gaussian model discussed 
in the text. 

i.e. a ~ y/~N for low-barrier ground states, well-described by GP theory and a ~ 1 in the 
fragmented limit. 

In Fig [TJ we show the how the value of (aja^) depends on a for a few decades of total 
particle number N. In Fig[7Ji), the failure of pu to stay at the correct asymptotic value of 
N/2 is due to the tails of the Gaussians exceeding the physical boundaries of the integration, 
a G {0, iV} and is an artifact of using this approximation for small particles numbers. 

Each of the graphs shows similar behavior as the system goes towards the fragmented 
limit. In order to see this in Fig [8] all four of the preceding graphs of pu has been scaled to 
unity and shown in the low-cr limit. 
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Since a two-mode model is only an appropriate description of a double-well condensate 
when the entanglement is less than unity, this analysis indicates that a large particle number 
BEC, starting in ground state of a harmonic trap, must be very precisely controlled in order 
to be a true two-mode system. Specifically, a ~ 2 has the physical meaning of a system which 
is well described by approximate Fock-states with an uncertainty in particle number of only 
a few particles since 2a ~ AiV. Condensates in splittin g ex periments performed to date, 
however, are typically on the order of ~ 10 5 particles j?! p\ js]. making it unlikely that 
such systems are in a regime where true two-mode physics is applicable for any significant 
amount of time. 

We have identified a universal feature of the splitting process for identical bosons, based 
solely on the assumption of a unimodal distribution of particles around a mean value, which 
indicates that independent of the total particle number or trap geometry, two-mode physics 
is only applicable when the uncertainty, AiV is on the order of a single particle. One of 
the primary implications of the iV-independence of this feature is that reliably scaling up 
experiments such as the lattice SF-MI experiments llj to much larger than unit filling of 
each well would require extremely fine control over the splitting potential. Further, figure 
[6] shows for an specific large- N system that the tunneling rate corresponding to two-mode 
physics is very small. 



VI. CONCLUSIONS 



We have presented a tractable method to describe the splitting of a large number of 
indistinguishable bosons in an external trap based on using localized non-linear functions 
and a novel mapping criterion to the low-barrier state instead of the typical approach of using 
linear combinations of non-localized functions. This method captures all of the features of 
the splitting of a BEC as a phase transition from a one-mode theory whose dynamical 
law is the GP equation, through a correlated two-mode theory and finally to a fragmented 
condensate defined by two independent complex order parameters, giving that system a 
U (1)@U (1) symmetry. By exploring the region which connects the one-mode and fragmented 
limits numerically and with an analytic model which captures the relevant behavior of the 
splitting process in Fock space, we are able to conclude that the region in which a true 
correlated two-mode model is applicable corresponds to non-vanishing fluctuations around 
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a well-defined number state of no more than one or two particles - more than this and full 
coherence is established between particles in the two wells and a single, de-localized, order 
parameter <£>op> whose time-evolution is dictated by the Gross-Pitaevskii equation, becomes 
a sufficient description of the state. 

We would like to thank David J. Masiello and Kaspar Sakmann for helpful discussions 
about multi-configurational methods and Nathan Kutz for sharing his expertise on numerical 
methods. We are also indebted to an anonymous referee of another draft for providing a set 
of detailed criteria for the clear presentation of a variationally-derived theory. We gratefully 
acknowledge funding from NSF grant PHY 07-03278. 

Appendix A: Effective parameters for 3D system 

We discuss a method to self-consistently include trap and mean-field physics from the 
condensate's extent in the two transverse directions when we are interested in the spatial 
profile (and ultimately dynamics) along only one principal splitting axis. This Appendix 
deals with producing effective lower- dimension equations of motion for NLSE/GP like equa- 
tions and multi-configurational systems of equations and is not concerned with systems in 
which the transverse confinement is on the order of the healing length or scattering length, 
at which point it has been shown [37] that a modified form of s-wave scattering takes place 
and an effective a s needs to be used. One can verify that none of the systems under consid- 
eration in this paper are in this limit, however, correct accounting of transverse degrees of 
freedom have an important effect on the results given above. 

For the sake of simplicity, we present this analysis by creating an effective ID equation 
for the x-direction from the full 3D GP equation, however it is exactly analogous to the 
treatment of the coupled mode equations above (I15p . which have all the same essential 
mathematical features of the GP equation. In three dimensions, the GP equation may be 
written: 



We start by assuming a separable solution of the GP equation ip(x,y, z) = <f>(x)ip±(y, z) 
with each factor independently normalized J dx\<fi(x)\ 2 = f dydz\ip±(y, z)\ 2 = 1, the time- 
independent 




(Al) 
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In the commonly-discussed experiment where a barrier is raised along the x-direction and 
the transverse directions are harmonically trapped, this becomes, under our assumption of 
separability: 

E<j>(x)Mv> z ) = ( V e*t(x) + \^ 2 y y 2 + \u>y - ^V 2 + g(N - l)|0(x)| 2 |^(y, z)\ 2 ^j 0(x)^±(y, z) 

(A2) 

Clearly, even under the assumption of a separable solution, the non-linear term does not 
allow the above equation to reduce to a set of two or three separable differential equations. 

While full 3D stationary and time-dependent solutions of the GP equation are regularly 
found {2 1 we want an effective ID equation with an eye towards performing tractable 
dynamical calculations at large particle numbers. The solution we propose is to assume 
a form for the transverse solution tpj_(y,z) which minimizes the energy of an appropriate 
energy functional given an initial solution of 4>{x), we then use a single value i/j±, such as 
the peak density, in the non-linear term of the - now separable - complex-time version of 
the GP equation which generates solutions of 4>(x). The new solution of <p(x) will give an 
updated energy functional for ip±(y,z) and the system of equations, eventually including 
the coupled mode and Fock-space equations for the two-mode theory presented in the body 
of this paper, can be iterated to convergence. 

Taking the Gaussian ansatz for the transverse profile 

1 =4 -* 2 

4>± (y, z) = e 2a v e^F (A3) 



V' 



We left-multiply flA2j) by ip(x,y,z) and integrate over all space, utilizing the fact that the 
x and transverse factors of ip are independently normalized in order to construct the full 
energy functional 

E = e x + e ± + g(N - l)r x r ± (A4) 

with the definitions e x = J dx<fi*(x)(V ex t(x) — ^jL)4>{x) 
^ = fdx<Kv,z)(±u}y* + - + &))Mv,*) 
and = fdx4>*(x)4>*(x)4>(x)(t)(x), 
T± = I dydzip*(y, z)^j*(y, z)<f>(y, z)(j)(y, z) 

The ansatz must minimize the energy E 1 - = e x + g(N — l)r x T or, explicitly 

E± = w (sj + r) + \ { « + ^ + 9(iV " »£r. (A5) 
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Where <j)(x) satisfies the effective ID equation 

M*) = (^ t (x)-^^ + ^(iV-l)^ 2 |^)| 2 )^) (A6) 

The two solutions can then be iterated to convergence through the parameters {i/j±, T 1 }, 
and provided that the dynamics in the transverse directions are negligible, yield effective 
ID time-dependent equations. 

Analogous to the case of the GP equation discussed above, the trap and particle data 
from the transverse directions can be accounted for in the two-mode theory discussed above. 
Summarizing the results from those calculations here ejj = ej- + e- 1 , = ej fc if j ^ k, 
Tjfcim = rj fcZm r- L . Further, in the nonlinear terms of the mode equations (|T5|) . a prefactor of 
tjj± 2 is included in both the diagonal and off-diagonal interaction strength. 

This procedure may be further generalized in a straight-forward manner to the case when 
dynamics in m dimensions are of interest, while the condensate is constrained to D > m 
dimensions. 
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